/* ptinpoly.c - point in polygon inside/outside code.

   by Eric Haines, 3D/Eye Inc, erich@eye.com

   This code contains the following algorithms:
    crossings - count the crossing made by a ray from the test point
    crossings-multiply - as above, but avoids a division; often a bit faster
    angle summation - sum the angle formed by point and vertex pairs
    weiler angle summation - sum the angles using quad movements
    half-plane testing - test triangle fan using half-space planes
    barycentric coordinates - test triangle fan w/barycentric coords
    spackman barycentric - preprocessed barycentric coordinates
    trapezoid testing - bin sorting algorithm
    grid testing - grid imposed on polygon
    exterior test - for convex polygons, check exterior of polygon
    inclusion test - for convex polygons, use binary search for edge.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ptinpoly.h"

#define X    0
#define Y    1

#ifndef TRUE
#define TRUE    1
#define FALSE   0
#endif

#ifndef HUGE
#define HUGE    1.797693134862315e+308
#endif

#ifndef M_PI
#define M_PI    3.14159265358979323846
#endif

/* test if a & b are within epsilon.  Favors cases where a < b */
#define Near(a,b,eps)    ( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) )

#define MALLOC_CHECK( a )    if ( !(a) ) {                   \
                    fprintf( stderr, "out of memory\n" );    \
                    exit(1);                                 \
                }


/* ======= Crossings algorithm ============================================ */

/* Shoot a test ray along +X axis.  The strategy, from MacMartin, is to
 * compare vertex Y values to the testing point's Y and quickly discard
 * edges which are entirely to one side of the test ray.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.    WINDING and CONVEX can be
 * defined for this test.
 */
int CrossingsTest(double pgon[][2], int numverts, double point[2])
{
#ifdef    WINDING
    int crossings;
#endif
    int j, yflag0, yflag1, inside_flag, xflag0;
    double ty, tx, * vtx0, * vtx1;
#ifdef    CONVEX
    int line_flag;
#endif

    tx = point[X];
    ty = point[Y];

    vtx0 = pgon[numverts - 1];
    /* get test bit for above/below X axis */
    yflag0 = (vtx0[Y] >= ty);
    vtx1 = pgon[0];

#ifdef    WINDING
    crossings = 0;
#else
    inside_flag = 0;
#endif
#ifdef    CONVEX
    line_flag = 0;
#endif
    for (j = numverts + 1; --j; ) {

        yflag1 = (vtx1[Y] >= ty);
        /* check if edge's endpoints straddle (are on opposite sides) of X axis
         * (i.e. the Y's differ); if so, +X ray could intersect this edge.
         */
        if (yflag0 != yflag1) {
            xflag0 = (vtx0[X] >= tx);
            /* check if endpoints are on same side of the Y axis (i.e. X's
             * are the same); if so, it's easy to test if edge hits or misses.
             */
            if (xflag0 == (vtx1[X] >= tx)) {

                /* if edge's X values both right of the point, must hit */
#ifdef    WINDING
                if (xflag0) crossings += (yflag0 ? -1 : 1);
#else
                if (xflag0) inside_flag = !inside_flag;
#endif
            }
            else {
                /* compute intersection of polygon edge with +X ray, note
                 * if >= point's X; if so, the ray hits it.
                 */
                if ((vtx1[X] - (vtx1[Y] - ty) *
                    (vtx0[X] - vtx1[X]) / (vtx0[Y] - vtx1[Y])) >= tx) {
#ifdef    WINDING
                    crossings += (yflag0 ? -1 : 1);
#else
                    inside_flag = !inside_flag;
#endif
                }
            }
#ifdef    CONVEX
            /* if this is second edge hit, then done testing */
            if (line_flag) goto Exit;

            /* note that one edge has been hit by the ray's line */
            line_flag = TRUE;
#endif
        }

        /* move to next pair of vertices, retaining info as possible */
        yflag0 = yflag1;
        vtx0 = vtx1;
        vtx1 += 2;
    }
#ifdef    CONVEX
    Exit : ;
#endif
#ifdef    WINDING
    /* test if crossings is not zero */
    inside_flag = (crossings != 0);
#endif

    return(inside_flag);
}

/* ======= Angle summation algorithm ======================================= */

/* Sum angles made by (vtxN to test point to vtxN+1), check if in proper
 * range to be inside.    VERY SLOW, included for tutorial reasons (i.e.
 * to show why one should never use this algorithm).
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 */
int AngleTest(double pgon[][2], int numverts, double point[2])
{
    double* vtx0, *vtx1, angle, len, vec0[2], vec1[2], vec_dot;
    int j;
    int inside_flag;

    /* sum the angles and see if answer mod 2*PI > PI */
    vtx0 = pgon[numverts - 1];
    vec0[X] = vtx0[X] - point[X];
    vec0[Y] = vtx0[Y] - point[Y];
    len = sqrt(vec0[X] * vec0[X] + vec0[Y] * vec0[Y]);
    if (len <= 0.0) {
        /* point and vertex coincide */
        return(1);
    }
    vec0[X] /= len;
    vec0[Y] /= len;

    angle = 0.0;
    for (j = 0; j < numverts; j++) {
        vtx1 = pgon[j];
        vec1[X] = vtx1[X] - point[X];
        vec1[Y] = vtx1[Y] - point[Y];
        len = sqrt(vec1[X] * vec1[X] + vec1[Y] * vec1[Y]);
        if (len <= 0.0) {
            /* point and vertex coincide */
            return(1);
        }
        vec1[X] /= len;
        vec1[Y] /= len;

        /* check if vec1 is to "left" or "right" of vec0 */
        vec_dot = vec0[X] * vec1[X] + vec0[Y] * vec1[Y];
        if (vec_dot < -1.0) {
            /* point is on edge, so always add 180 degrees.  always
             * adding is not necessarily the right answer for
             * concave polygons and subtractive triangles.
             */
            angle += M_PI;
        }
        else if (vec_dot < 1.0) {
            if (vec0[X] * vec1[Y] - vec1[X] * vec0[Y] >= 0.0) {
                /* add angle due to dot product of vectors */
                angle += acos(vec_dot);
            }
            else {
                /* subtract angle due to dot product of vectors */
                angle -= acos(vec_dot);
            }
        } /* if vec_dot >= 1.0, angle does not change */

        /* get to next point */
        vtx0 = vtx1;
        vec0[X] = vec1[X];
        vec0[Y] = vec1[Y];
    }
    /* test if between PI and 3*PI, 5*PI and 7*PI, etc */
    inside_flag = fmod(fabs(angle) + M_PI, 4.0 * M_PI) > 2.0 * M_PI;

    return(inside_flag);
}

/* ======= Weiler algorithm ============================================ */

/* Track quadrant movements using Weiler's algorithm (elsewhere in Graphics
 * Gems IV).  Algorithm has been optimized for testing purposes, but the
 * crossings algorithm is still faster. Included to provide timings.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 *
 * WINDING can be defined for this test.
 */

#define QUADRANT( vtx, x, y )    \
    (((vtx)[X]>(x)) ? ( ((vtx)[Y]>(y)) ? 0:3 ) : ( ((vtx)[Y]>(y)) ? 1:2 ))

#define X_INTERCEPT( v0, v1, y )    \
    ( (((v1)[X]-(v0)[X])/((v1)[Y]-(v0)[Y])) * ((y)-(v0)[Y]) + (v0)[X] )

int WeilerTest(double pgon[][2], int numverts, double point[2])
{
    int angle, qd, next_qd, delta, j;
    double ty, tx, *vtx0, *vtx1;
    int inside_flag;

    tx = point[X];
    ty = point[Y];

    vtx0 = pgon[numverts - 1];
    qd = QUADRANT(vtx0, tx, ty);
    angle = 0;

    vtx1 = pgon[0];

    for (j = numverts + 1; --j; ) {
        /* calculate quadrant and delta from last quadrant */
        next_qd = QUADRANT(vtx1, tx, ty);
        delta = next_qd - qd;

        /* adjust delta and add it to total angle sum */
        switch (delta) {
        case 0:    /* do nothing */
            break;
        case -1:
        case 3:
            angle--;
            qd = next_qd;
            break;

        case 1:
        case -3:
            angle++;
            qd = next_qd;
            break;

        case 2:
        case -2:
            if (X_INTERCEPT(vtx0, vtx1, ty) > tx) {
                angle -= delta;
            }
            else {
                angle += delta;
            }
            qd = next_qd;
            break;
        }

        /* increment for next step */
        vtx0 = vtx1;
        vtx1 += 2;
    }

#ifdef    WINDING
    /* simple windings test:  if angle != 0, then point is inside */
    inside_flag = (angle != 0);
#else
    /* Jordan test:  if angle is +-4, 12, 20, etc then point is inside */
    inside_flag = ((angle / 4) & 0x1);
#endif

    return (inside_flag);
}

#undef    QUADRANT
#undef    Y_INTERCEPT

/* ======= Triangle half-plane algorithm ================================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point is inside of the three half-planes formed by the triangle's edges.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to a plane set array.
 * Call testing procedure with a pointer to this array, _numverts_, and
 * test point _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to plane set array to free space.
 *
 * SORT and CONVEX can be defined for this test.
 */

 /* split polygons along set of x axes - call preprocess once */
pPlaneSet PlaneSetup(double pgon[][2], int numverts)
{
    int i, p1, p2;
    double tx, ty, vx0, vy0;
    pPlaneSet pps, pps_return;
#ifdef    SORT
    double len[3], len_temp;
    int j;
    PlaneSet    ps_temp;
#ifdef    CONVEX
    pPlaneSet    pps_new;
    pSizePlanePair p_size_pair;
#endif
#endif

    pps = pps_return =
        (pPlaneSet)malloc(3 * (numverts - 2) * sizeof(PlaneSet));
    MALLOC_CHECK(pps);
#ifdef    CONVEX
#ifdef    SORT
    p_size_pair =
        (pSizePlanePair)malloc((numverts - 2) * sizeof(SizePlanePair));
    MALLOC_CHECK(p_size_pair);
#endif
#endif

    vx0 = pgon[0][X];
    vy0 = pgon[0][Y];

    for (p1 = 1, p2 = 2; p2 < numverts; p1++, p2++) {
        pps->vx = vy0 - pgon[p1][Y];
        pps->vy = pgon[p1][X] - vx0;
        pps->c = pps->vx * vx0 + pps->vy * vy0;
#ifdef    SORT
        len[0] = pps->vx * pps->vx + pps->vy * pps->vy;
#ifdef    CONVEX
#ifdef    HYBRID
        pps->ext_flag = (p1 == 1);
#endif
        /* Sort triangles by areas, so compute (twice) the area here */
        p_size_pair[p1 - 1].pps = pps;
        p_size_pair[p1 - 1].size =
            (pgon[0][X] * pgon[p1][Y]) +
            (pgon[p1][X] * pgon[p2][Y]) +
            (pgon[p2][X] * pgon[0][Y]) -
            (pgon[p1][X] * pgon[0][Y]) -
            (pgon[p2][X] * pgon[p1][Y]) -
            (pgon[0][X] * pgon[p2][Y]);
#endif
#endif
        pps++;
        pps->vx = pgon[p1][Y] - pgon[p2][Y];
        pps->vy = pgon[p2][X] - pgon[p1][X];
        pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y];
#ifdef    SORT
        len[1] = pps->vx * pps->vx + pps->vy * pps->vy;
#endif
#ifdef    CONVEX
#ifdef    HYBRID
        pps->ext_flag = TRUE;
#endif
#endif
        pps++;
        pps->vx = pgon[p2][Y] - vy0;
        pps->vy = vx0 - pgon[p2][X];
        pps->c = pps->vx * pgon[p2][X] + pps->vy * pgon[p2][Y];
#ifdef    SORT
        len[2] = pps->vx * pps->vx + pps->vy * pps->vy;
#endif
#ifdef    CONVEX
#ifdef    HYBRID
        pps->ext_flag = (p2 == numverts - 1);
#endif
#endif

        /* find an average point which must be inside of the triangle */
        tx = (vx0 + pgon[p1][X] + pgon[p2][X]) / 3.0;
        ty = (vy0 + pgon[p1][Y] + pgon[p2][Y]) / 3.0;

        /* check sense and reverse if test point is not thought to be inside
         * first triangle
         */
        if (pps->vx * tx + pps->vy * ty >= pps->c) {
            /* back up to start of plane set */
            pps -= 2;
            /* point is thought to be outside, so reverse sense of edge
             * normals so that it is correctly considered inside.
             */
            for (i = 0; i < 3; i++) {
                pps->vx = -pps->vx;
                pps->vy = -pps->vy;
                pps->c = -pps->c;
                pps++;
            }
        }
        else {
            pps++;
        }

#ifdef    SORT
        /* sort the planes based on the edge lengths */
        pps -= 3;
        for (i = 0; i < 2; i++) {
            for (j = i + 1; j < 3; j++) {
                if (len[i] < len[j]) {
                    ps_temp = pps[i];
                    pps[i] = pps[j];
                    pps[j] = ps_temp;
                    len_temp = len[i];
                    len[i] = len[j];
                    len[j] = len_temp;
                }
            }
        }
        pps += 3;
#endif
    }

#ifdef    CONVEX
#ifdef    SORT
    /* sort the triangles based on their areas */
    qsort(p_size_pair, numverts - 2,
        sizeof(SizePlanePair), CompareSizePlanePairs);

    /* make the plane sets match the sorted order */
    for (i = 0, pps = pps_return
        ; i < numverts - 2
        ; i++) {

        pps_new = p_size_pair[i].pps;
        for (j = 0; j < 3; j++, pps++, pps_new++) {
            ps_temp = *pps;
            *pps = *pps_new;
            *pps_new = ps_temp;
        }
    }
    free(p_size_pair);
#endif
#endif

    return(pps_return);
}

#ifdef    CONVEX
#ifdef    SORT
int CompareSizePlanePairs(const void *p0, const void* p1)
{
    pSizePlanePair p_sp0 = (pSizePlanePair)p0;
    pSizePlanePair p_sp1 = (pSizePlanePair)p1;
    if (p_sp0->size == p_sp1->size) {
        return(0);
    }
    else {
        return(p_sp0->size > p_sp1->size ? -1 : 1);
    }
}
#endif
#endif


/* check point for inside of three "planes" formed by triangle edges */
int PlaneTest(pPlaneSet p_plane_set, int numverts, double point[2])
{
    pPlaneSet ps;
    int p2;
#ifndef CONVEX
    int inside_flag;
#endif
    double tx, ty;

    tx = point[X];
    ty = point[Y];

#ifndef CONVEX
    inside_flag = 0;
#endif

    for (ps = p_plane_set, p2 = numverts - 1; --p2; ) {

        if (ps->vx * tx + ps->vy * ty < ps->c) {
            ps++;
            if (ps->vx * tx + ps->vy * ty < ps->c) {
                ps++;
                /* note: we make the third edge have a slightly different
                 * equality condition, since this third edge is in fact
                 * the next triangle's first edge.  Not fool-proof, but
                 * it doesn't hurt (better would be to keep track of the
                 * triangle's area sign so we would know which kind of
                 * triangle this is).  Note that edge sorting nullifies
                 * this special inequality, too.
                 */
                if (ps->vx * tx + ps->vy * ty <= ps->c) {
                    /* point is inside polygon */
#ifdef CONVEX
                    return(1);
#else
                    inside_flag = !inside_flag;
#endif
                }
#ifdef    CONVEX
#ifdef    HYBRID
                /* check if outside exterior edge */
                else if (ps->ext_flag) return(0);
#endif
#endif
                ps++;
            }
            else {
#ifdef    CONVEX
#ifdef    HYBRID
                /* check if outside exterior edge */
                if (ps->ext_flag) return(0);
#endif
#endif
                /* get past last two plane tests */
                ps += 2;
            }
        }
        else {
#ifdef    CONVEX
#ifdef    HYBRID
            /* check if outside exterior edge */
            if (ps->ext_flag) return(0);
#endif
#endif
            /* get past all three plane tests */
            ps += 3;
        }
    }

#ifdef CONVEX
    /* for convex, if we make it to here, all triangles were missed */
    return(0);
#else
    return(inside_flag);
#endif
}

void PlaneCleanup(pPlaneSet p_plane_set)
{
    free(p_plane_set);
}

/* ======= Barycentric algorithm ========================================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point has barycentric coordinates which are inside the triangle.
 * Similar to Badouel's code in Graphics Gems, with a little more efficient
 * coding.
 *
 * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
 * _point_, returns 1 if inside, 0 if outside.
 */
int BarycentricTest(double pgon[][2], int numverts, double point[2])
{
    double* pg1, *pg2, *pgend;
    double tx, ty, u0, u1, u2, v0, v1, vx0, vy0, alpha, beta, denom;
    int inside_flag;

    tx = point[X];
    ty = point[Y];
    vx0 = pgon[0][X];
    vy0 = pgon[0][Y];
    u0 = tx - vx0;
    v0 = ty - vy0;

    inside_flag = 0;
    pgend = pgon[numverts - 1];
    for (pg1 = pgon[1], pg2 = pgon[2]; pg1 != pgend; pg1 += 2, pg2 += 2) {

        u1 = pg1[X] - vx0;
        if (u1 == 0.0) {

            /* 0 and 1 vertices have same X value */

            /* zero area test - can be removed for convex testing */
            u2 = pg2[X] - vx0;
            if ((u2 == 0.0) ||

                /* compute beta and check bounds */
                /* we use "<= 0.0" so that points on the shared interior
                 * edge will (generally) be inside only one polygon.
                 */
                ((beta = u0 / u2) <= 0.0) ||
                (beta > 1.0) ||

                /* zero area test - remove for convex testing */
                ((v1 = pg1[Y] - vy0) == 0.0) ||

                /* compute alpha and check bounds */
                ((alpha = (v0 - beta *
                (pg2[Y] - vy0)) / v1) < 0.0)) {

                /* whew! missed! */
                goto NextTri;
            }

        }
        else {
            /* 0 and 1 vertices have different X value */

            /* compute denom and check for zero area triangle - check
             * is not needed for convex polygon testing
             */
            u2 = pg2[X] - vx0;
            v1 = pg1[Y] - vy0;
            denom = (pg2[Y] - vy0) * u1 - u2 * v1;
            if ((denom == 0.0) ||

                /* compute beta and check bounds */
                /* we use "<= 0.0" so that points on the shared interior
                 * edge will (generally) be inside only one polygon.
                 */
                ((beta = (v0 * u1 - u0 * v1) / denom) <= 0.0) ||
                (beta > 1.0) ||

                /* compute alpha & check bounds */
                ((alpha = (u0 - beta * u2) / u1) < 0.0)) {

                /* whew! missed! */
                goto NextTri;
            }
        }

        /* check gamma */
        if (alpha + beta <= 1.0) {
            /* survived */
            inside_flag = !inside_flag;
        }

    NextTri:;
    }
    return(inside_flag);
}

/* ======= Barycentric precompute (Spackman) algorithm ===================== */

/* Split the polygon into a fan of triangles and for each triangle test if
 * the point has barycentric coordinates which are inside the triangle.
 * Use Spackman's normalization method to precompute various parameters.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to the array of the parameters records and the
 * number of parameter records created.
 * Call testing procedure with the first vertex in the polygon _pgon[0]_,
 * a pointer to this array, the number of parameter records, and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to parameter record array to free space.
 *
 * SORT can be defined for this test.
 * (CONVEX could be added: see PlaneSetup and PlaneTest for method)
 */
pSpackmanSet SpackmanSetup(double pgon[][2], int numverts, int* p_numrec)
{
    int p1, p2, degen;
    double denom, u1, v1, *pv[3];
    pSpackmanSet pss, pss_return;
#ifdef    SORT
    double u[2], v[2], len[2], *pv_temp;
#endif

    pss = pss_return =
        (pSpackmanSet)malloc((numverts - 2) * sizeof(SpackmanSet));
    MALLOC_CHECK(pss);

    degen = 0;

    for (p1 = 1, p2 = 2; p2 < numverts; p1++, p2++) {

        pv[0] = pgon[0];
        pv[1] = pgon[p1];
        pv[2] = pgon[p2];

#ifdef    SORT
        /* Note that sorting can cause a mismatch of alpha/beta inequality
         * tests.  In other words, test points on an interior line between
         * test triangles will often then be wrong.
         */
        u[0] = pv[1][X] - pv[0][X];
        u[1] = pv[2][X] - pv[0][X];
        v[0] = pv[1][Y] - pv[0][Y];
        v[1] = pv[2][Y] - pv[0][Y];
        len[0] = u[0] * u[0] + v[0] * v[0];
        len[1] = u[1] * u[1] + v[1] * v[1];

        /* compare two edges touching anchor point and put longest first */
        /* we don't sort all three edges because the anchor point and
         * values computed from it gets used for all triangles in the fan.
         */
        if (len[0] < len[1]) {
            pv_temp = pv[1];
            pv[1] = pv[2];
            pv[2] = pv_temp;
        }
#endif

        u1 = pv[1][X] - pv[0][X];
        pss->u2 = pv[2][X] - pv[0][X];
        v1 = pv[1][Y] - pv[0][Y];
        pss->v2 = pv[2][Y] - pv[0][Y];
        pss->u1_nonzero = !(u1 == 0.0);
        if (pss->u1_nonzero) {
            /* not zero, so compute inverse */
            pss->inv_u1 = 1.0 / u1;
            denom = pss->v2 * u1 - pss->u2 * v1;
            if (denom == 0.0) {
                /* degenerate triangle, ignore it */
                degen++;
                goto Skip;
            }
            else {
                pss->u1p = u1 / denom;
                pss->v1p = v1 / denom;
            }
        }
        else {
            if (pss->u2 == 0.0) {
                /* degenerate triangle, ignore it */
                degen++;
                goto Skip;
            }
            else {
                /* not zero, so compute inverse */
                pss->inv_u2 = 1.0 / pss->u2;
                if (v1 == 0.0) {
                    /* degenerate triangle, ignore it */
                    degen++;
                    goto Skip;
                }
                else {
                    pss->inv_v1 = 1.0 / v1;
                }
            }
        }

        pss++;
    Skip:;
    }

    /* number of Spackman records */
    *p_numrec = numverts - degen - 2;
    if (degen) {
        pss = pss_return =
            (pSpackmanSet)realloc(pss_return,
            (numverts - 2 - degen) * sizeof(SpackmanSet));
    }

    return(pss_return);
}

/* barycentric, a la Badouel's Graphics Gems article and Spackman's normalization precompute */
int SpackmanTest(double anchor[2], pSpackmanSet p_spackman_set, int numrec, double point[2])
{
    pSpackmanSet pss;
    int inside_flag;
    int nr;
    double tx, ty, vx0, vy0, u0, v0, alpha, beta;

    tx = point[X];
    ty = point[Y];
    /* note that we really need only the first vertex of the polygon,
     * so do not really need to keep the whole polygon around.
     */
    vx0 = anchor[X];
    vy0 = anchor[Y];
    u0 = tx - vx0;
    v0 = ty - vy0;

    inside_flag = 0;

    for (pss = p_spackman_set, nr = numrec + 1; --nr; pss++) {

        if (pss->u1_nonzero) {
            /* 0 and 2 vertices have different X value */

            /* compute beta and check bounds */
            /* we use "<= 0.0" so that points on the shared interior edge
             * will (generally) be inside only one polygon.
             */
            beta = (v0 * pss->u1p - u0 * pss->v1p);
            if ((beta <= 0.0) || (beta > 1.0) ||

                /* compute alpha & check bounds */
                ((alpha = (u0 - beta * pss->u2) * pss->inv_u1)
                    < 0.0)) {

                /* whew! missed! */
                goto NextTri;
            }
        }
        else {
            /* 0 and 2 vertices have same X value */

            /* compute beta and check bounds */
            /* we use "<= 0.0" so that points on the shared interior edge
             * will (generally) be inside only one polygon.
             */
            beta = u0 * pss->inv_u2;
            if ((beta <= 0.0) || (beta >= 1.0) ||

                /* compute alpha and check bounds */
                ((alpha = (v0 - beta * pss->v2) * pss->inv_v1)
                    < 0.0)) {

                /* whew! missed! */
                goto NextTri;
            }
        }

        /* check gamma */
        if (alpha + beta <= 1.0) {
            /* survived */
            inside_flag = !inside_flag;
        }

    NextTri:;
    }

    return(inside_flag);
}

void SpackmanCleanup(pSpackmanSet p_spackman_set)
{
    free(p_spackman_set);
}

/* ======= Trapezoid (bin) algorithm ======================================= */

/* Split polygons along set of y bins and sorts the edge fragments.  Testing
 * is done against these fragments.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices, the
 * number of bins desired _bins_, and a pointer to a trapezoid structure
 * _p_trap_set_.
 * Call testing procedure with 2D polygon _pgon_ with _numverts_ number of
 * vertices, _p_trap_set_ pointer to trapezoid structure, and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to trapezoid structure to free space.
 */
void TrapezoidSetup(double pgon[][2], int numverts, int bins, pTrapezoidSet p_trap_set)
{
    double* vtx0, *vtx1, *vtxa, *vtxb, slope;
    int i, j, bin_tot[TOT_BINS], ba, bb, id, full_cross, count;
    double fba, fbb, vx0, vx1, dy, vy0;

    p_trap_set->bins = bins;
    p_trap_set->trapz = (pTrapezoid)malloc(p_trap_set->bins *
        sizeof(Trapezoid));
    MALLOC_CHECK(p_trap_set->trapz);

    p_trap_set->minx =
        p_trap_set->maxx = pgon[0][X];
    p_trap_set->miny =
        p_trap_set->maxy = pgon[0][Y];

    for (i = 1; i < numverts; i++) {
        if (p_trap_set->minx > (vx0 = pgon[i][X])) {
            p_trap_set->minx = vx0;
        }
        else if (p_trap_set->maxx < vx0) {
            p_trap_set->maxx = vx0;
        }

        if (p_trap_set->miny > (vy0 = pgon[i][Y])) {
            p_trap_set->miny = vy0;
        }
        else if (p_trap_set->maxy < vy0) {
            p_trap_set->maxy = vy0;
        }
    }

    /* add a little to the bounds to ensure everything falls inside area */
    p_trap_set->miny -= EPSILON * (p_trap_set->maxy - p_trap_set->miny);
    p_trap_set->maxy += EPSILON * (p_trap_set->maxy - p_trap_set->miny);

    p_trap_set->ydelta =
        (p_trap_set->maxy - p_trap_set->miny) / (double)p_trap_set->bins;
    p_trap_set->inv_ydelta = 1.0 / p_trap_set->ydelta;

    /* find how many locations to allocate for each bin */
    for (i = 0; i < p_trap_set->bins; i++) {
        bin_tot[i] = 0;
    }

    vtx0 = pgon[numverts - 1];
    for (i = 0; i < numverts; i++) {
        vtx1 = pgon[i];

        /* skip if Y's identical (edge has no effect) */
        if (vtx0[Y] != vtx1[Y]) {

            if (vtx0[Y] < vtx1[Y]) {
                vtxa = vtx0;
                vtxb = vtx1;
            }
            else {
                vtxa = vtx1;
                vtxb = vtx0;
            }
            ba = (int)((vtxa[Y] - p_trap_set->miny) * p_trap_set->inv_ydelta);
            fbb = (vtxb[Y] - p_trap_set->miny) * p_trap_set->inv_ydelta;
            bb = (int)fbb;
            /* if high vertex ends on a boundary, don't go into next boundary */
            if (fbb - (double)bb == 0.0) {
                bb--;
            }

            /* mark the bins with this edge */
            for (j = ba; j <= bb; j++) {
                bin_tot[j]++;
            }
        }

        vtx0 = vtx1;
    }

    /* allocate the bin contents and fill in some basics */
    for (i = 0; i < p_trap_set->bins; i++) {
        p_trap_set->trapz[i].edge_set =
            (pEdge*)malloc(bin_tot[i] * sizeof(pEdge));
        MALLOC_CHECK(p_trap_set->trapz[i].edge_set);
        for (j = 0; j < bin_tot[i]; j++) {
            p_trap_set->trapz[i].edge_set[j] =
                (pEdge)malloc(sizeof(Edge));
            MALLOC_CHECK(p_trap_set->trapz[i].edge_set[j]);
        }

        /* start these off at some awful values; refined below */
        p_trap_set->trapz[i].minx = p_trap_set->maxx;
        p_trap_set->trapz[i].maxx = p_trap_set->minx;
        p_trap_set->trapz[i].count = 0;
    }

    /* now go through list yet again, putting edges in bins */
    vtx0 = pgon[numverts - 1];
    id = numverts - 1;
    for (i = 0; i < numverts; i++) {
        vtx1 = pgon[i];

        /* we can skip edge if Y's are equal */
        if (vtx0[Y] != vtx1[Y]) {
            if (vtx0[Y] < vtx1[Y]) {
                vtxa = vtx0;
                vtxb = vtx1;
            }
            else {
                vtxa = vtx1;
                vtxb = vtx0;
            }
            fba = (vtxa[Y] - p_trap_set->miny) * p_trap_set->inv_ydelta;
            ba = (int)fba;
            fbb = (vtxb[Y] - p_trap_set->miny) * p_trap_set->inv_ydelta;
            bb = (int)fbb;
            /* if high vertex ends on a boundary, don't go into it */
            if (fbb == (double)bb) {
                bb--;
            }

            vx0 = vtxa[X];
            dy = vtxa[Y] - vtxb[Y];
            slope = p_trap_set->ydelta * (vtxa[X] - vtxb[X]) / dy;

            /* set vx1 in case loop is not entered */
            vx1 = vx0;
            full_cross = 0;

            for (j = ba; j < bb; j++, vx0 = vx1) {
                /* could increment vx1, but for greater accuracy recompute it */
                vx1 = vtxa[X] + ((double)(j + 1) - fba) * slope;

                count = p_trap_set->trapz[j].count++;
                p_trap_set->trapz[j].edge_set[count]->id = id;
                p_trap_set->trapz[j].edge_set[count]->full_cross = full_cross;
                TrapBound(j, count, vx0, vx1, p_trap_set);
                full_cross = 1;
            }

            /* at last bin - fill as above, but with vx1 = vtxb[X] */
            vx0 = vx1;
            vx1 = vtxb[X];
            count = p_trap_set->trapz[bb].count++;
            p_trap_set->trapz[bb].edge_set[count]->id = id;
            /* the last bin is never a full crossing */
            p_trap_set->trapz[bb].edge_set[count]->full_cross = 0;
            TrapBound(bb, count, vx0, vx1, p_trap_set);
        }

        vtx0 = vtx1;
        id = i;
    }

    /* finally, sort the bins' contents by minx */
    for (i = 0; i < p_trap_set->bins; i++) {
        qsort(p_trap_set->trapz[i].edge_set, p_trap_set->trapz[i].count,
            sizeof(pEdge), CompareEdges);
    }
}

void TrapBound(int j, int count, double vx0, double vx1, pTrapezoidSet p_trap_set)
{
    double xt;

    if (vx0 > vx1) {
        xt = vx0;
        vx0 = vx1;
        vx1 = xt;
    }

    if (p_trap_set->trapz[j].minx > vx0) {
        p_trap_set->trapz[j].minx = vx0;
    }
    if (p_trap_set->trapz[j].maxx < vx1) {
        p_trap_set->trapz[j].maxx = vx1;
    }
    p_trap_set->trapz[j].edge_set[count]->minx = vx0;
    p_trap_set->trapz[j].edge_set[count]->maxx = vx1;
}

/* used by qsort to sort */
int CompareEdges(const void* uu, const void* vv)
{
    pEdge* u = (pEdge*)uu;
    pEdge* v = (pEdge*)vv;

    if ((*u)->minx == (*v)->minx) {
        return(0);
    }
    else {
        return((*u)->minx < (*v)->minx ? -1 : 1);
    }
}

int TrapezoidTest(double pgon[][2], int  numverts, pTrapezoidSet    p_trap_set, double point[2])
{
    int j, b, count, id;
    double tx, ty, *vtx0, *vtx1;
    pEdge* pp_bin;
    pTrapezoid p_trap;
    int inside_flag;

    inside_flag = 0;

    /* first, is point inside bounding rectangle? */
    if ((ty = point[Y]) < p_trap_set->miny ||
        ty >= p_trap_set->maxy ||
        (tx = point[X]) < p_trap_set->minx ||
        tx >= p_trap_set->maxx) {

        /* outside of box */
        return(0);
    }

    /* what bin are we in? */
    b = (int)((ty - p_trap_set->miny) * p_trap_set->inv_ydelta);

    /* find if we're inside this bin's bounds */
    if (tx < (p_trap = &p_trap_set->trapz[b])->minx ||
        tx > p_trap->maxx) {

        /* outside of box */
        return(0);
    }

    /* now search bin for crossings */
    pp_bin = p_trap->edge_set;
    count = p_trap->count;
    for (j = 0; j < count; j++, pp_bin++) {
        if (tx < (*pp_bin)->minx) {

            /* all remaining edges are to right of point, so test them */
            do {
                if ((*pp_bin)->full_cross) {
                    inside_flag = !inside_flag;
                }
                else {
                    id = (*pp_bin)->id;
                    if ((ty <= pgon[id][Y]) !=
                        (ty <= pgon[(id + 1) % numverts][Y])) {

                        /* point crosses edge in Y, so must cross */
                        inside_flag = !inside_flag;
                    }
                }
                pp_bin++;
            } while (++j < count);
            goto Exit;

        }
        else if (tx < (*pp_bin)->maxx) {
            /* edge is overlapping point in X, check it */
            id = (*pp_bin)->id;
            vtx0 = pgon[id];
            vtx1 = pgon[(id + 1) % numverts];

            if ((*pp_bin)->full_cross ||
                (ty <= vtx0[Y]) != (ty <= vtx1[Y])) {

                /* edge crosses in Y, so have to do full crossings test */
                if ((vtx0[X] -
                    (vtx0[Y] - ty) *
                    (vtx1[X] - vtx0[X]) / (vtx1[Y] - vtx0[Y])) >= tx) {
                    inside_flag = !inside_flag;
                }
            }

        } /* else edge is to left of point, ignore it */
    }

Exit:
    return(inside_flag);
}

void TrapezoidCleanup(pTrapezoidSet p_trap_set)
{
    int i, j, count;

    for (i = 0; i < p_trap_set->bins; i++) {
        /* all of these should have bin sets, but check just in case */
        if (p_trap_set->trapz[i].edge_set) {
            count = p_trap_set->trapz[i].count;
            for (j = 0; j < count; j++) {
                if (p_trap_set->trapz[i].edge_set[j]) {
                    free(p_trap_set->trapz[i].edge_set[j]);
                }
            }
            free(p_trap_set->trapz[i].edge_set);
        }
    }
    free(p_trap_set->trapz);
}

/* ======= Grid algorithm ================================================= */

/* Impose a grid upon the polygon and test only the local edges against the
 * point.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * grid resolution _resolution_ and a pointer to a grid structure _p_gs_.
 * Call testing procedure with a pointer to this array and test point _point_,
 * returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to grid structure to free space.
 */

 /* Strategy for setup:
  *   Get bounds of polygon, allocate grid.
  *   "Walk" each edge of the polygon and note which edges have been crossed
  *     and what cells are entered (points on a grid edge are always considered
  *     to be above that edge).    Keep a record of the edges overlapping a cell.
  *     For cells with edges, determine if any cell border has no edges passing
  *     through it and so can be used for shooting a test ray.
  *     Keep track of the parity of the x (horizontal) grid cell borders for
  *     use in determining whether the grid corners are inside or outside.
  */
void GridSetup(double pgon[][2], int numverts, int resolution, pGridSet p_gs)
{
    double *vtx0, *vtx1, *vtxa, *vtxb, *p_gl;
    int i, j, gc_clear_flags;
    double vx0, vx1, vy0, vy1, gxdiff, gydiff, eps;
    pGridCell p_gc, p_ngc;
    double xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty;
    double tgcx, tgcy;
    int gcx, gcy, sign_x;
    int y_flag, io_state;

    p_gs->xres = p_gs->yres = resolution;
    p_gs->tot_cells = p_gs->xres * p_gs->yres;
    p_gs->glx = (double*)malloc((p_gs->xres + 1) * sizeof(double));
    MALLOC_CHECK(p_gs->glx);
    p_gs->gly = (double*)malloc((p_gs->yres + 1) * sizeof(double));
    MALLOC_CHECK(p_gs->gly);
    p_gs->gc = (pGridCell)malloc(p_gs->tot_cells * sizeof(GridCell));
    MALLOC_CHECK(p_gs->gc);

    p_gs->minx =
        p_gs->maxx = pgon[0][X];
    p_gs->miny =
        p_gs->maxy = pgon[0][Y];

    /* find bounds of polygon */
    for (i = 1; i < numverts; i++) {
        vx0 = pgon[i][X];
        if (p_gs->minx > vx0) {
            p_gs->minx = vx0;
        }
        else if (p_gs->maxx < vx0) {
            p_gs->maxx = vx0;
        }

        vy0 = pgon[i][Y];
        if (p_gs->miny > vy0) {
            p_gs->miny = vy0;
        }
        else if (p_gs->maxy < vy0) {
            p_gs->maxy = vy0;
        }
    }

    /* add a little to the bounds to ensure everything falls inside area */
    gxdiff = p_gs->maxx - p_gs->minx;
    gydiff = p_gs->maxy - p_gs->miny;
    p_gs->minx -= EPSILON * gxdiff;
    p_gs->maxx += EPSILON * gxdiff;
    p_gs->miny -= EPSILON * gydiff;
    p_gs->maxy += EPSILON * gydiff;

    /* avoid roundoff problems near corners by not getting too close to them */
    eps = 1e-9 * (gxdiff + gydiff);

    /* use the new bounds to compute cell widths */
TryAgain:
    p_gs->xdelta =
        (p_gs->maxx - p_gs->minx) / (double)p_gs->xres;
    p_gs->inv_xdelta = 1.0 / p_gs->xdelta;

    p_gs->ydelta =
        (p_gs->maxy - p_gs->miny) / (double)p_gs->yres;
    p_gs->inv_ydelta = 1.0 / p_gs->ydelta;

    for (i = 0, p_gl = p_gs->glx; i < p_gs->xres; i++) {
        *p_gl++ = p_gs->minx + i * p_gs->xdelta;
    }
    /* make last grid corner precisely correct */
    *p_gl = p_gs->maxx;

    for (i = 0, p_gl = p_gs->gly; i < p_gs->yres; i++) {
        *p_gl++ = p_gs->miny + i * p_gs->ydelta;
    }
    *p_gl = p_gs->maxy;

    for (i = 0, p_gc = p_gs->gc; i < p_gs->tot_cells; i++, p_gc++) {
        p_gc->tot_edges = 0;
        p_gc->gc_flags = 0x0;
        p_gc->gr = NULL;
    }

    /* loop through edges and insert into grid structure */
    vtx0 = pgon[numverts - 1];
    for (i = 0; i < numverts; i++) {
        vtx1 = pgon[i];

        if (vtx0[Y] < vtx1[Y]) {
            vtxa = vtx0;
            vtxb = vtx1;
        }
        else {
            vtxa = vtx1;
            vtxb = vtx0;
        }

        /* Set x variable for the direction of the ray */
        xdiff = vtxb[X] - vtxa[X];
        ydiff = vtxb[Y] - vtxa[Y];
        tmax = sqrt(xdiff * xdiff + ydiff * ydiff);

        /* if edge is of 0 length, ignore it (useless edge) */
        if (tmax == 0.0) goto NextEdge;

        xdir = xdiff / tmax;
        ydir = ydiff / tmax;

        gcx = (int)((vtxa[X] - p_gs->minx) * p_gs->inv_xdelta);
        gcy = (int)((vtxa[Y] - p_gs->miny) * p_gs->inv_ydelta);

        /* get information about slopes of edge, etc */
        if (vtxa[X] == vtxb[X]) {
            sign_x = 0;
            tx = HUGE;
            tgcx = 0;
        }
        else {
            inv_x = tmax / xdiff;
            tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X];
            if (vtxa[X] < vtxb[X]) {
                sign_x = 1;
                tx += p_gs->xdelta;
                tgcx = p_gs->xdelta * inv_x;
            }
            else {
                sign_x = -1;
                tgcx = -p_gs->xdelta * inv_x;
            }
            tx *= inv_x;
        }

        if (vtxa[Y] == vtxb[Y]) {
            ty = HUGE;
            tgcy = 0;
        }
        else {
            inv_y = tmax / ydiff;
            ty = (p_gs->ydelta * (double)(gcy + 1) + p_gs->miny - vtxa[Y])
                * inv_y;
            tgcy = p_gs->ydelta * inv_y;
        }

        p_gc = &p_gs->gc[gcy * p_gs->xres + gcx];

        vx0 = vtxa[X];
        vy0 = vtxa[Y];

        t_near = 0.0;

        do {
            /* choose the next boundary, but don't move yet */
            if (tx <= ty) {
                gcx += sign_x;

                ty -= tx;
                t_near += tx;
                tx = tgcx;

                /* note which edge is hit when leaving this cell */
                if (t_near < tmax) {
                    if (sign_x > 0) {
                        p_gc->gc_flags |= GC_R_EDGE_HIT;
                        vx1 = p_gs->glx[gcx];
                    }
                    else {
                        p_gc->gc_flags |= GC_L_EDGE_HIT;
                        vx1 = p_gs->glx[gcx + 1];
                    }

                    /* get new location */
                    vy1 = t_near * ydir + vtxa[Y];
                }
                else {
                    /* end of edge, so get exact value */
                    vx1 = vtxb[X];
                    vy1 = vtxb[Y];
                }

                y_flag = FALSE;

            }
            else {

                gcy++;

                tx -= ty;
                t_near += ty;
                ty = tgcy;

                /* note top edge is hit when leaving this cell */
                if (t_near < tmax) {
                    p_gc->gc_flags |= GC_T_EDGE_HIT;
                    /* this toggles the parity bit */
                    p_gc->gc_flags ^= GC_T_EDGE_PARITY;

                    /* get new location */
                    vx1 = t_near * xdir + vtxa[X];
                    vy1 = p_gs->gly[gcy];
                }
                else {
                    /* end of edge, so get exact value */
                    vx1 = vtxb[X];
                    vy1 = vtxb[Y];
                }

                y_flag = TRUE;
            }

            /* check for corner crossing, then mark the cell we're in */
            if (!AddGridRecAlloc(p_gc, vx0, vy0, vx1, vy1, eps)) {
                /* warning, danger - we have just crossed a corner.
                 * There are all kinds of topological messiness we could
                 * do to get around this case, but they're a headache.
                 * The simplest recovery is just to change the extents a bit
                 * and redo the meshing, so that hopefully no edges will
                 * perfectly cross a corner.  Since it's a preprocess, we
                 * don't care too much about the time to do it.
                 */

                 /* clean out all grid records */
                for (i = 0, p_gc = p_gs->gc
                    ; i < p_gs->tot_cells
                    ; i++, p_gc++) {

                    if (p_gc->gr) {
                        free(p_gc->gr);
                    }
                }

                /* make the bounding box ever so slightly larger, hopefully
                 * changing the alignment of the corners.
                 */
                p_gs->minx -= EPSILON * gxdiff * 0.24;
                p_gs->miny -= EPSILON * gydiff * 0.10;

                /* yes, it's the dreaded goto - run in fear for your lives! */
                goto TryAgain;
            }

            if (t_near < tmax) {
                /* note how we're entering the next cell */
                /* TBD: could be done faster by incrementing index in the
                 * incrementing code, above */
                p_gc = &p_gs->gc[gcy * p_gs->xres + gcx];

                if (y_flag) {
                    p_gc->gc_flags |= GC_B_EDGE_HIT;
                    /* this toggles the parity bit */
                    p_gc->gc_flags ^= GC_B_EDGE_PARITY;
                }
                else {
                    p_gc->gc_flags |=
                        (sign_x > 0) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT;
                }
            }

            vx0 = vx1;
            vy0 = vy1;
        }
        /* have we gone further than the end of the edge? */
        while (t_near < tmax);

    NextEdge:
        vtx0 = vtx1;
    }

    /* the grid is all set up, now set up the inside/outside value of each
     * corner.
     */
    p_gc = p_gs->gc;
    p_ngc = &p_gs->gc[p_gs->xres];

    /* we know the bottom and top rows are all outside, so no flag is set */
    for (i = 1; i < p_gs->yres; i++) {
        /* start outside */
        io_state = 0x0;

        for (j = 0; j < p_gs->xres; j++) {

            if (io_state) {
                /* change cell left corners to inside */
                p_gc->gc_flags |= GC_TL_IN;
                p_ngc->gc_flags |= GC_BL_IN;
            }

            if (p_gc->gc_flags & GC_T_EDGE_PARITY) {
                io_state = !io_state;
            }

            if (io_state) {
                /* change cell right corners to inside */
                p_gc->gc_flags |= GC_TR_IN;
                p_ngc->gc_flags |= GC_BR_IN;
            }

            p_gc++;
            p_ngc++;
        }
    }

    p_gc = p_gs->gc;
    for (i = 0; i < p_gs->tot_cells; i++) {

        /* reverse parity of edge clear (1==edge clear) */
        gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR;
        if (gc_clear_flags & GC_L_EDGE_CLEAR) {
            p_gc->gc_flags |= GC_AIM_L;
        }
        else
            if (gc_clear_flags & GC_B_EDGE_CLEAR) {
                p_gc->gc_flags |= GC_AIM_B;
            }
            else
                if (gc_clear_flags & GC_R_EDGE_CLEAR) {
                    p_gc->gc_flags |= GC_AIM_R;
                }
                else
                    if (gc_clear_flags & GC_T_EDGE_CLEAR) {
                        p_gc->gc_flags |= GC_AIM_T;
                    }
                    else {
                        /* all edges have something on them, do full test */
                        p_gc->gc_flags |= GC_AIM_C;
                    }
        p_gc++;
    }
}

int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps)
{
    pGridRec p_gr;
    double slope, inv_slope;

    if (Near(ya, yb, eps)) {
        if (Near(xa, xb, eps)) {
            /* edge is 0 length, so get rid of it */
            return(FALSE);
        }
        else {
            /* horizontal line */
            slope = HUGE;
            inv_slope = 0.0;
        }
    }
    else {
        if (Near(xa, xb, eps)) {
            /* vertical line */
            slope = 0.0;
            inv_slope = HUGE;
        }
        else {
            slope = (xb - xa) / (yb - ya);
            inv_slope = (yb - ya) / (xb - xa);
        }
    }

    p_gc->tot_edges++;
    if (p_gc->tot_edges <= 1) {
        p_gc->gr = (pGridRec)malloc(sizeof(GridRec));
    }
    else {
        p_gc->gr = (pGridRec)realloc(p_gc->gr,
            p_gc->tot_edges * sizeof(GridRec));
    }
    MALLOC_CHECK(p_gc->gr);
    p_gr = &p_gc->gr[p_gc->tot_edges - 1];

    p_gr->slope = slope;
    p_gr->inv_slope = inv_slope;

    p_gr->xa = xa;
    p_gr->ya = ya;
    if (xa <= xb) {
        p_gr->minx = xa;
        p_gr->maxx = xb;
    }
    else {
        p_gr->minx = xb;
        p_gr->maxx = xa;
    }
    if (ya <= yb) {
        p_gr->miny = ya;
        p_gr->maxy = yb;
    }
    else {
        p_gr->miny = yb;
        p_gr->maxy = ya;
    }

    /* P2 - P1 */
    p_gr->ax = xb - xa;
    p_gr->ay = yb - ya;

    return(TRUE);
}

/* Test point against grid and edges in the cell (if any).  Algorithm:
 *    Check bounding box; if outside then return.
 *    Check cell point is inside; if simple inside or outside then return.
 *    Find which edge or corner is considered to be the best for testing and
 *      send a test ray towards it, counting the crossings.  Add in the
 *      state of the edge or corner the ray went to and so determine the
 *      state of the point (inside or outside).
 */
int GridTest(pGridSet p_gs, double point[2])
{
    int j, count, init_flag;
    pGridCell p_gc;
    pGridRec p_gr;
    double tx, ty, xcell, ycell, bx, by, cx, cy, cornerx, cornery;
    double alpha, beta, denom;
    unsigned short gc_flags;
    int inside_flag = FALSE;

    /* first, is point inside bounding rectangle? */
    if ((ty = point[Y]) < p_gs->miny ||
        ty >= p_gs->maxy ||
        (tx = point[X]) < p_gs->minx ||
        tx >= p_gs->maxx) {

        /* outside of box */
        //inside_flag = FALSE;
    }
    else {

        /* what cell are we in? */
        ycell = (ty - p_gs->miny) * p_gs->inv_ydelta;
        xcell = (tx - p_gs->minx) * p_gs->inv_xdelta;
        p_gc = &p_gs->gc[((int)ycell) * p_gs->xres + (int)xcell];

        /* is cell simple? */
        count = p_gc->tot_edges;
        if (count) {
            /* no, so find an edge which is free. */
            gc_flags = p_gc->gc_flags;
            p_gr = p_gc->gr;
            switch (gc_flags & GC_AIM) {
            case GC_AIM_L:
                /* left edge is clear, shoot X- ray */
                /* note - this next statement requires that GC_BL_IN is 1 */
                inside_flag = gc_flags & GC_BL_IN;
                for (j = count + 1; --j; p_gr++) {
                    /* test if y is between edges */
                    if (ty >= p_gr->miny && ty < p_gr->maxy) {
                        if (tx > p_gr->maxx) {
                            inside_flag = !inside_flag;
                        }
                        else if (tx > p_gr->minx) {
                            /* full computation */
                            if ((p_gr->xa -
                                (p_gr->ya - ty) * p_gr->slope) < tx) {
                                inside_flag = !inside_flag;
                            }
                        }
                    }
                }
                break;

            case GC_AIM_B:
                /* bottom edge is clear, shoot Y+ ray */
                /* note - this next statement requires that GC_BL_IN is 1 */
                inside_flag = gc_flags & GC_BL_IN;
                for (j = count + 1; --j; p_gr++) {
                    /* test if x is between edges */
                    if (tx >= p_gr->minx && tx < p_gr->maxx) {
                        if (ty > p_gr->maxy) {
                            inside_flag = !inside_flag;
                        }
                        else if (ty > p_gr->miny) {
                            /* full computation */
                            if ((p_gr->ya - (p_gr->xa - tx) *
                                p_gr->inv_slope) < ty) {
                                inside_flag = !inside_flag;
                            }
                        }
                    }
                }
                break;

            case GC_AIM_R:
                /* right edge is clear, shoot X+ ray */
                inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0;

                /* TBD: Note, we could have sorted the edges to be tested
                 * by miny or somesuch, and so be able to cut testing
                 * short when the list's miny > point.y .
                 */
                for (j = count + 1; --j; p_gr++) {
                    /* test if y is between edges */
                    if (ty >= p_gr->miny && ty < p_gr->maxy) {
                        if (tx <= p_gr->minx) {
                            inside_flag = !inside_flag;
                        }
                        else if (tx <= p_gr->maxx) {
                            /* full computation */
                            if ((p_gr->xa -
                                (p_gr->ya - ty) * p_gr->slope) >= tx) {
                                inside_flag = !inside_flag;
                            }
                        }
                    }
                }
                break;

            case GC_AIM_T:
                /* top edge is clear, shoot Y+ ray */
                inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0;
                for (j = count + 1; --j; p_gr++) {
                    /* test if x is between edges */
                    if (tx >= p_gr->minx && tx < p_gr->maxx) {
                        if (ty <= p_gr->miny) {
                            inside_flag = !inside_flag;
                        }
                        else if (ty <= p_gr->maxy) {
                            /* full computation */
                            if ((p_gr->ya - (p_gr->xa - tx) *
                                p_gr->inv_slope) >= ty) {
                                inside_flag = !inside_flag;
                            }
                        }
                    }
                }
                break;

            case GC_AIM_C:
                /* no edge is clear, bite the bullet and test
                 * against the bottom left corner.
                 * We use Franklin Antonio's algorithm (Graphics Gems III).
                 */
                 /* TBD: Faster yet might be to test against the closest
                  * corner to the cell location, but our hope is that we
                  * rarely need to do this testing at all.
                  */
                inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN);
                init_flag = TRUE;

                /* get lower left corner coordinate */
                cornerx = p_gs->glx[(int)xcell];
                cornery = p_gs->gly[(int)ycell];
                for (j = count + 1; --j; p_gr++) {

                    /* quick out test: if test point is
                     * less than minx & miny, edge cannot overlap.
                     */
                    if (tx >= p_gr->minx && ty >= p_gr->miny) {

                        /* quick test failed, now check if test point and
                         * corner are on different sides of edge.
                         */
                        if (init_flag) {
                            /* Compute these at most once for test */
                            /* P3 - P4 */
                            bx = tx - cornerx;
                            by = ty - cornery;
                            init_flag = FALSE;
                        }
                        /* you may get a warning about bx and by not being initialized, but they are, above */
                        denom = p_gr->ay * bx - p_gr->ax * by;
                        if (denom != 0.0) {
                            /* lines are not collinear, so continue */
                            /* P1 - P3 */
                            cx = p_gr->xa - tx;
                            cy = p_gr->ya - ty;
                            alpha = by * cx - bx * cy;
                            if (denom > 0.0) {
                                if (alpha < 0.0 || alpha >= denom) {
                                    /* test edge not hit */
                                    goto NextEdge;
                                }
                                beta = p_gr->ax * cy - p_gr->ay * cx;
                                if (beta < 0.0 || beta >= denom) {
                                    /* polygon edge not hit */
                                    goto NextEdge;
                                }
                            }
                            else {
                                if (alpha > 0.0 || alpha <= denom) {
                                    /* test edge not hit */
                                    goto NextEdge;
                                }
                                beta = p_gr->ax * cy - p_gr->ay * cx;
                                if (beta > 0.0 || beta <= denom) {
                                    /* polygon edge not hit */
                                    goto NextEdge;
                                }
                            }
                            inside_flag = !inside_flag;
                        }

                    }
                NextEdge:;
                }
                break;
            }
        }
        else {
            /* simple cell, so if lower left corner is in,
             * then cell is inside.
             */
            inside_flag = p_gc->gc_flags & GC_BL_IN;
        }
    }

    return(inside_flag);
}

void GridCleanup(pGridSet p_gs)
{
    int i;
    pGridCell p_gc;

    for (i = 0, p_gc = p_gs->gc
        ; i < p_gs->tot_cells
        ; i++, p_gc++) {

        if (p_gc->gr) {
            free(p_gc->gr);
        }
    }
    free(p_gs->glx);
    free(p_gs->gly);
    free(p_gs->gc);
}

/* ======= Exterior (convex only) algorithm =============================== */

/* Test the edges of the convex polygon against the point.  If the point is
 * outside any edge, the point is outside the polygon.
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to a plane set array.
 * Call testing procedure with a pointer to this array, _numverts_, and
 * test point _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to plane set array to free space.
 *
 * RANDOM can be defined for this test.
 * CONVEX must be defined for this test; it is not usable for general polygons.
 */

#ifdef    CONVEX
 /* make exterior plane set */
pPlaneSet ExteriorSetup(double pgon[][2], int numverts)
{
    int p1, p2, flip_edge;
    pPlaneSet pps, pps_return;
#ifdef    RANDOM
    int i, ind;
    PlaneSet ps_temp;
#endif

    pps = pps_return =
        (pPlaneSet)malloc(numverts * sizeof(PlaneSet));
    MALLOC_CHECK(pps);

    /* take cross product of vertex to find handedness */
    flip_edge = (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y]) >
        (pgon[0][Y] - pgon[1][Y])* (pgon[1][X] - pgon[2][X]);

    /* Generate half-plane boundary equations now for faster testing later.
     * vx & vy are the edge's normal, c is the offset from the origin.
     */
    for (p1 = numverts - 1, p2 = 0; p2 < numverts; p1 = p2, p2++, pps++) {
        pps->vx = pgon[p1][Y] - pgon[p2][Y];
        pps->vy = pgon[p2][X] - pgon[p1][X];
        pps->c = pps->vx * pgon[p1][X] + pps->vy * pgon[p1][Y];

        /* check sense and reverse plane edge if need be */
        if (flip_edge) {
            pps->vx = -pps->vx;
            pps->vy = -pps->vy;
            pps->c = -pps->c;
        }
    }

#ifdef    RANDOM
    /* Randomize the order of the edges to improve chance of early out */
    /* There are better orders, but the default order is the worst */
    for (i = 0, pps = pps_return
        ; i < numverts
        ; i++) {

        ind = (int)(RAN01() * numverts);
        if ((ind < 0) || (ind >= numverts)) {
            fprintf(stderr,
                "Yikes, the random number generator is returning values\n");
            fprintf(stderr,
                "outside the range [0.0,1.0), so please fix the code!\n");
            ind = 0;
        }

        /* swap edges */
        ps_temp = *pps;
        *pps = pps_return[ind];
        pps_return[ind] = ps_temp;
    }
#endif

    return(pps_return);
}

/* Check point for outside of all planes */
/* note that we don't need "pgon", since it's been processed into
 * its corresponding PlaneSet.
 */
int ExteriorTest(pPlaneSet p_ext_set, int numverts, double point[2])
{
    PlaneSet* pps;
    int p0;
    double tx, ty;

    tx = point[X];
    ty = point[Y];

    for (p0 = numverts + 1, pps = p_ext_set; --p0; pps++) {

        /* test if the point is outside this edge */
        if (pps->vx * tx + pps->vy * ty > pps->c) {
            return(0);
        }
    }
    /* if we make it to here, we were inside all edges */
    return(1);
}

void ExteriorCleanup(pPlaneSet p_ext_set)
{
    free(p_ext_set);
}
#endif

/* ======= Inclusion (convex only) algorithm ============================== */

/* Create an efficiency structure (see Preparata) for the convex polygon which
 * allows binary searching to find which edge to test the point against.  This
 * algorithm is O(log n).
 *
 * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices,
 * which returns a pointer to an inclusion anchor structure.
 * Call testing procedure with a pointer to this structure and test point
 * _point_, returns 1 if inside, 0 if outside.
 * Call cleanup with pointer to inclusion anchor structure to free space.
 *
 * CONVEX must be defined for this test; it is not usable for general polygons.
 */

#ifdef    CONVEX
 /* make inclusion wedge set */
pInclusionAnchor InclusionSetup(double pgon[][2], int numverts)
{
    int pc, p1, p2, flip_edge;
    double ax, ay, qx, qy, wx, wy, len;
    pInclusionAnchor pia;
    pInclusionSet pis;

    /* double the first edge to avoid needing modulo during test search */
    pia = (pInclusionAnchor)malloc(sizeof(InclusionAnchor));
    MALLOC_CHECK(pia);
    pis = pia->pis =
        (pInclusionSet)malloc((numverts + 1) * sizeof(InclusionSet));
    MALLOC_CHECK(pis);

    pia->hi_start = numverts - 1;

    /* get average point to make wedges from */
    qx = qy = 0.0;
    for (p2 = 0; p2 < numverts; p2++) {
        qx += pgon[p2][X];
        qy += pgon[p2][Y];
    }
    pia->qx = qx /= (double)numverts;
    pia->qy = qy /= (double)numverts;

    /* take cross product of vertex to find handedness */
    pia->flip_edge = flip_edge =
        (pgon[0][X] - pgon[1][X]) * (pgon[1][Y] - pgon[2][Y]) >
        (pgon[0][Y] - pgon[1][Y])* (pgon[1][X] - pgon[2][X]);


    ax = pgon[0][X] - qx;
    ay = pgon[0][Y] - qy;
    len = sqrt(ax * ax + ay * ay);
    if (len == 0.0) {
        fprintf(stderr, "sorry, polygon for inclusion test is defective\n");
        exit(1);
    }
    pia->ax = ax /= len;
    pia->ay = ay /= len;

    /* loop through edges, and double last edge */
    for (pc = p1 = 0, p2 = 1
        ; pc <= numverts
        ; pc++, p1 = p2, p2 = (++p2) % numverts, pis++) {

        /* wedge border */
        wx = pgon[p1][X] - qx;
        wy = pgon[p1][Y] - qy;
        len = sqrt(wx * wx + wy * wy);
        wx /= len;
        wy /= len;

        /* cosine of angle from anchor border to wedge border */
        pis->dot = ax * wx + ay * wy;
        /* sign from cross product */
        if ((ax * wy > ay* wx) == flip_edge) {
            pis->dot = -2.0 - pis->dot;
        }

        /* edge */
        pis->ex = pgon[p1][Y] - pgon[p2][Y];
        pis->ey = pgon[p2][X] - pgon[p1][X];
        pis->ec = pis->ex * pgon[p1][X] + pis->ey * pgon[p1][Y];

        /* check sense and reverse plane eqns if need be */
        if (flip_edge) {
            pis->ex = -pis->ex;
            pis->ey = -pis->ey;
            pis->ec = -pis->ec;
        }
    }
    /* set first angle a little > 1.0 and last < -3.0 just to be safe. */
    pia->pis[0].dot = -3.001;
    pia->pis[numverts].dot = 1.001;

    return(pia);
}

/* Find wedge point is in by binary search, then test wedge */
int InclusionTest(pInclusionAnchor pia, double point[2])
{
    double tx, ty, len, dot;
    int inside_flag, lo, hi, ind;
    pInclusionSet pis;

    tx = point[X] - pia->qx;
    ty = point[Y] - pia->qy;
    len = sqrt(tx * tx + ty * ty);
    /* check if point is exactly at anchor point (which is inside polygon) */
    if (len == 0.0) return(1);
    tx /= len;
    ty /= len;

    /* get dot product for searching */
    dot = pia->ax * tx + pia->ay * ty;
    if ((pia->ax * ty > pia->ay* tx) == pia->flip_edge) {
        dot = -2.0 - dot;
    }

    /* binary search through angle list and find matching angle pair */
    lo = 0;
    hi = pia->hi_start;
    while (lo <= hi) {
        ind = (lo + hi) / 2;
        if (dot < pia->pis[ind].dot) {
            hi = ind - 1;
        }
        else if (dot > pia->pis[ind + 1].dot) {
            lo = ind + 1;
        }
        else {
            goto Foundit;
        }
    }
    /* should never reach here, but just in case... */
    fprintf(stderr,
        "Hmmm, something weird happened - bad dot product %lg\n", dot);

Foundit:

    /* test if the point is outside the wedge's exterior edge */
    /* may get a warning that ind may not be initialized; it should be */
    pis = &pia->pis[ind];
    inside_flag = (pis->ex * point[X] + pis->ey * point[Y] <= pis->ec);

    return(inside_flag);
}

void InclusionCleanup(pInclusionAnchor p_inc_anchor)
{
    free(p_inc_anchor->pis);
    free(p_inc_anchor);
}
#endif


/* ======= Crossings Multiply algorithm =================================== */

/*
 * This version is usually somewhat faster than the original published in
 * Graphics Gems IV; by turning the division for testing the X axis crossing
 * into a tricky multiplication test this part of the test became faster,
 * which had the additional effect of making the test for "both to left or
 * both to right" a bit slower for triangles than simply computing the
 * intersection each time.  The main increase is in triangle testing speed,
 * which was about 15% faster; all other polygon complexities were pretty much
 * the same as before.  On machines where division is very expensive (not the
 * case on the HP 9000 series on which I tested) this test should be much
 * faster overall than the old code.  Your mileage may (in fact, will) vary,
 * depending on the machine and the test data, but in general I believe this
 * code is both shorter and faster.  This test was inspired by unpublished
 * Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
 * Related work by Samosky is in:
 *
 * Samosky, Joseph, "SectionView: A system for interactively specifying and
 * visualizing sections through three-dimensional medical image data",
 * M.S. Thesis, Department of Electrical Engineering and Computer Science,
 * Massachusetts Institute of Technology, 1993.
 *
 */

 /* Shoot a test ray along +X axis.  The strategy is to compare vertex Y values
  * to the testing point's Y and quickly discard edges which are entirely to one
  * side of the test ray.  Note that CONVEX and WINDING code can be added as
  * for the CrossingsTest() code; it is left out here for clarity.
  *
  * Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
  * _point_, returns 1 if inside, 0 if outside.
  */
int CrossingsMultiplyTest(double pgon[][2], int numverts, double point[2])
{
    int j, yflag0, yflag1, inside_flag;
    double ty, tx, *vtx0, *vtx1;

    tx = point[X];
    ty = point[Y];

    vtx0 = pgon[numverts - 1];
    /* get test bit for above/below X axis */
    yflag0 = (vtx0[Y] >= ty);
    vtx1 = pgon[0];

    inside_flag = 0;
    for (j = numverts + 1; --j; ) {

        yflag1 = (vtx1[Y] >= ty);
        /* Check if endpoints straddle (are on opposite sides) of X axis
         * (i.e. the Y's differ); if so, +X ray could intersect this edge.
         * The old test also checked whether the endpoints are both to the
         * right or to the left of the test point.  However, given the faster
         * intersection point computation used below, this test was found to
         * be a break-even proposition for most polygons and a loser for
         * triangles (where 50% or more of the edges which survive this test
         * will cross quadrants and so have to have the X intersection computed
         * anyway).  I credit Joseph Samosky with inspiring me to try dropping
         * the "both left or both right" part of my code.
         */
        if (yflag0 != yflag1) {
            /* Check intersection of pgon segment with +X ray.
             * Note if >= point's X; if so, the ray hits it.
             * The division operation is avoided for the ">=" test by checking
             * the sign of the first vertex wrto the test point; idea inspired
             * by Joseph Samosky's and Mark Haigh-Hutchinson's different
             * polygon inclusion tests.
             */
            if (((vtx1[Y] - ty) * (vtx0[X] - vtx1[X]) >=
                (vtx1[X] - tx) * (vtx0[Y] - vtx1[Y])) == yflag1) {
                inside_flag = !inside_flag;
            }
        }

        /* Move to the next pair of vertices, retaining info as possible. */
        yflag0 = yflag1;
        vtx0 = vtx1;
        vtx1 += 2;
    }

    return(inside_flag);
}
